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1. INTRODUCTION 

In the biological activated sludge process (ASP), microbial contaminants are used in influent 
wastewater to degrade organic matters and nutrients [1-3]. Due to the system’s dynamic uncertainties 
and highly nonlinearities in the microorganism’s growth kinetic model, control of activated sludge process 
is a challenge. Conventional controllers are designed based on mathematical model of ASP which requires 
the designer to spend lot of building the model from the basic laws of biomass mass balance. However, there 
are model free design methods available that can be utilized in designing control systems [4-6]. 
For conventional control, the differential equations derived are then transformed into a standard transfer 
function or state space and the control laws are then applied [7-10]. However, the principal control problems 
of ASP as found in literature can be briefly illustrated by the variability of the kinetic parameter [11-12], 
time-varying influent conditions, nonlinearities [13], delays, lack of accurate sensors [14], and the limited 
availability of the online measurements [15]. Therefore, adaptive control principles are the preferred 
design option [16-18]. 

The main issue in the design of predictive controller for ASP was how to implement an accurate 
internal predictive model which best describes the behavior of the biochemical process involved. 
In [19] the internal predictive model of the plant implemented neural networks for controller design. 
However, this method suffers from the fact that it requires much attention on the training data collection 
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resulting in many approximations in this stage of the design, and huge validation experiments have to be 
done to accurately identify the model and it’s usually fails to track the variability of the kinetic parameters. 
In [20] the authors have built the internal model for prediction based on linearized state space model of 
the bioprocess which is definitely undesirable in cases of variability in the effluents or any changes in 
operating conditions. Subspace Identification techniques developed in nineties offers an attractive 
and effective solution for the issue of prediction model and has been deployed successfully in design of 
predictive controllers for multiple-input multiple-output processes [21]. 

The employment of the robust computational tools such as QR-decomposition has recently opened 
a wide area of research to develop a single step implementation for predictive controllers where only 
input-output data are processed to directly calculate controller parameters or namely “Model-free” control 
methods [22]. The results of [22-24] are some of the successful reports for this class of novel algorithms. 
However, they have not introduced any closed-loop stability analysis for the sake of practical implementation 
of the proposed controller. Most adaptive predictive controllers reported so far are only applied to linear 
systems or indirectly implemented. Therefore, the stability involved can be easily examined and proved for 
several cases of infinite prediction horizon or when the open loop system can be accessed and is stable. 
However, when the system is nonlinear or when finite horizons is applied, only few approaches have been 
published and were limited to indirect schemes [25-26]. In this paper, a data-driven predictive controller for 
activated sludge process is applied which is also consider stability of the proposed controller. 


2. ACTIVATED SLUDGE MODEL DEVELOPMENT 

The biological process is modeled from component-balance considerations and the obtained model 
is most often a highly complex and high-order non-linear system given by the following set of nonlinear 
differential equations [27]: 


dX(t) 


O Z WOX) -DOC + HX) + DOX-C) (1) 
BO Z -EOX — D+ SC) + DO)Sin (2) 
AO = — SHO x(t) — DO + NCU) + Kia(Cs — C(t) + DOHCin (3) 
SO = DO + DX(H) — DOG + NX (4) 


where X(t), S(t), X,(t), C(t) are the state variables representing the biomass, the substrate, the recycled 
biomass and dissolved oxygen concentrations, respectively. Meanwhile, D(t) is the dilution rate, and r and B 
represent the ratio of recycled flow to influent flow and ratio of waste flow to influent flow, Sj, and Cin 
respectively, which corresponds to the substrate and dissolved oxygen concentrations in the influent stream. 
In addition, Y is the cell mass; Ky is a constant term, C, is the maximum dissolved oxygen concentration 
and Ka represents the oxygen mass transfer coefficient. 


Kya = aW, a>0 (5) 


where W is the airflow rate. The kinetic model of the specific cell mass production rate p(t) is given as [24]. 


= s® cW 
L(t) = Umax ks KCO (6) 








where [pax is the maximum specific growth rate, K, is the so-called affinity constant, expressing 
the dependency of the degradation rate on the concentration of pollutant S, and K, is the saturation constant. 
Notice that all the above equations are only used here in this work to generate input-output (I/O) data 
and they will not be used in any stage of the controller design. For simulations purpose, the process 
parameters, kinetic parameters, and initial conditions are adopted, as shown in Tables 1-3, respectively. 








able 1. Process parameters able 2. Kinetic parameters able 3. Initial conditions 
Table 1. P p t Table 2. Kinetic p t Table 3. Initial dit 
Parameter Value Parameter Value Parameter Value 

Y 0.65 Y 0.65 X(0) 215 mgl” 

r 0.6 T 0.6 S(0) 35 mel" 

B 0.2 B 0.2 C(O) 6 mgl” 

a 0.018 m° a 0.018 m° Xy(0) 400 mg!" 

Ko 0.5 Ko 0.5 Sin 200 mel! 

C; 10 mgl! Cy 10 mgl! Cin 0.5 mgl’ 

Sin 200 mgt! Sin 200 mgt! 

Cin 0.5 mgl’ Cin 0.5 mgl’ 
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3. SUBSPACE BASED PREDICTOR IDENTIFICATION 
The most commonly used technique in system identification is used in this paper, given by (7): 


X1 = Ax, + Bu, + Ke; and y, = Cx, + Duy + ef (7) 


where, x; € R”, is the system states with order n, u E€ Ru, is the input vector with N, number of inputs, 
Yt € RYY is the output vector with Ny number of outputs and, e, is a white noise disturbance (innovations) 
sequence with covariance Ze. The matrices A € R™™", B € R™*Nu, C € R™ Ny, D € RNu*Ny and K € RNy*Ny 
are the state, input, output, feed-forward and Kalman gain matrices of the system. The following (8) to (10) 
are defined as fundamental subspace I/O matrix equations. The subscripts p here denotes “past” and f 
denotes “future” matrices for its respective variables. 


= d 
Yp = TiXp + H@U, + HEE, (8) 
Ye = TiX_ + HoUy + HSE¢ (9) 
Xp = AX, + AGU, + ASE, (10) 


where the extended observability matrix I, the reversed extended controllability matrix for deterministic 
input Ad and for stochastic input A}, matrices H and HẸ are the lower triangular Toeplitz matrices for 
both deterministic and stochastic parts, respectively. Assuming that the open-loop I/O data u, and y, for 
k = {0,1,...,N — 1} are available for the identification, define the past input, U, and future input, Uç data 


matrices for (8)-(10) as follows: 


Uo Uy vt UN-2i 

Up # Upon #| 2 U2 UN (11) 
Ui-1 Uj UN-i 
Uj Ui+1 1 UyN-i 

Ur £ Uäijzi-1) Ê oe oe E het (12) 
Ugi-1 Uži ss UN-1 


where {Up, Us} € IRNuixN-2i_ The subscripts (k,|k,) indicates that the first column of the subspace matrix 
starts from the time step k, and ends by the time step kz. Similar notation can be written for past and future 
output matrices. In subspace identification literature, the following short-hand notation is commonly used: 


Wp = B (13) 


Ye = LwWp + LU; (14) 


Equation (14) is the subspace predictor equation where the parameters to be obtained contain system model 
information without completely identifying the parametric model of the system such as the conventional 
transfer function or state space representation of the system. In addition, L,,RiNy*i(Ny+Nu). is the subspace 
matrix that corresponds to the past inputs and outputs. L,R'Ny*!Nu: is the subspace matrix that corresponds to 
the future input data or the deterministic inputs, and Ly = Hê. The problem of obtaining the best linear 
predictor of the future outputs Y;, given the past inputs and outputs in the form of Wp and the future inputs Us 


can be written as a Frobenius norm minimization: 


(15) 








Ye Mw Lall] 








2 
MIN Ly Lu 
F 


The least square problem in (15) can be solved by the orthogonal projection and can be implemented by 


forming a QR-decomposition of the matrix A Uy yY | ‘as follows: 
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W] Ri 0 0] 
Us| =|R2 R22 0 JQ. (16) 
Yr R31 R32 R334 LQ3 











Thus, we can estimate optimal prediction Y; as described in (16). Ly, and L, can be used to formulate 
an estimation of future outputs as follows: 





Uj 
Vx ig Uk 
A : = Lw] Yk-i Abg : (17) 
Vkti-1 : Uk-i-1 
Yk-1 


Equation (17) is a streamlined subspace-based linear predictor equation which can be written as, 


Ve = LwWp + LyU¢ (18) 


4. DIRECT ADAPTIVE MODEL PREDICTIVE CONTROL 

In the direct adaptive model predictive control (DAMPC) or so-called data driven method, 
the controller parameters are estimated directly from the measurements. Therefore, the system model 
is implicitly included and no explicit model is formulated at any stage of the implementation. The adaptive 
mechanism implemented here adjusts the predictive controller parameters to the changes in the process 
dynamics in every sampling time. The measured input-output data is collected over a sliding window. It can 
be noted that the data window used to identify the controller parameters or state matrices can be expressed in 
terms of future inputs, us = [Ui-1 ‘** Uzi-2]" and measurement (past) inputs Up = [Yo +] uj-1]" 
and outputs yp = [Yo '] Yi-1]". Here, two prediction problems are given at current time instant i 
and i+1 as shown in Figure 2. The first prediction problem (t = i) represents the case for obtaining 
the optimal prediction of i future outputs f= [Yj - Fezail using the information given in 
the previously stated data window up, Yp and ur. The second prediction problem shows that the time instant 
slides from t = ito t= i+ 1, this differs from the data window (up, Yp and uf) which is now slides from left 
to right. At every time step, for the new available input-output data, the MPC parameters are updated online 
according to the updated predictor gain matrices and a new control action is computed. The cost function 
performance criterion used can be written in matrix form as, 


J = Or- rA" Q r- ro + Au, R Aun, (19) 
Iny Aut+1 
In Auts2 $ A $ zags 
where, rp =] .” |r > Aun, = : . By incorporating the integral action to the subspace predictive 
Iny Aui+H, 
model, (20) is derived: 
Şe = Fy Yt + Fy, LwAWp + Fy, by Aug, (20) 
In, = 0 
z ‘ H Inen : 
where, Fy =] i ™ : JandLyo =Ly | a |. Then, we can rewrite (19) as: 
y On xHe 
Ny e Iny 
J= Z Auk QAU, + Au} = (21) 


with, 0 = (Fu, La) TQ (Fp, Li) +R and E = (Fy, Lie) TQ (Fa,LwAw, + Fy, (Ye - ma) In this case, 


the model predictive law for the unconstrained case of the problem of (21) is derived, which can be simply 
solved by: 
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sy > 
Tagg = MA, +E = 0 
Aup, = —0718 


At time t, only Au(1) is implemented and the calculation is repeated at each time instant. Hence, the control 
law is as follows: 


Au(1) = —Kaw,AWp = Ke (Yt — rt+1) (22) 


with, Kawp = Un, ONyx(He-1) Ny] Kawp He and Ke = Un, ONux(He-)Nul Ken, Where Kawp.He and Ken, 


defined as, 
Kawptte = ((FnyLu’) 7Q (Fx, La’) + R) (Fuy Lu) TQ (Fn, Lw) 
Ken, = ((Fn, Li) 1O (Fn, L) + R) (Fu L) TO Fn, 


5. SIMULATION AND RESULTS 
The implemented control system is a two inputs two outputs system that are dilution rate, D and air 


flow rate, W for the inputs, and substrate, S and dissolved oxygen, C for the outputs. The simulation is run 
from a steady state operating point at outputs S=41.23 mg/l, C=6.11 mg/l and inputs D=0.081/hr and W=90 
m3/hr. The sampling period is chosen as Ts=1. The weighting matrices are chosen as Q = diag(1,10) 
and R = diag(10*,10~*). Meanwhile, the Prediction Horizon, H, and Control Horizon, H, were 35,and 5 
respectively. Figure 1 shows the evolution of the outputs and inputs. The setpoint given for the outputs are 
varied for an approximately 10% around the system’s steady state condition. The squared error performance 
index was chosen to evaluate the performance of the closed-loop system and it was found to be 445.2. 


Substrate Concentration Response 
44 
S —— DAMPC 
eee Ba pe S O aoa Na Re 1] 
es 0 
4 
1200 1300 1400 1500 1600 1700 1800 1900 2000 
Dissolved Oxygen Concentration Response 
65 
> —— DAMPC 
x Le ee ee = sf hd eal ee Ref | 
6 
1200 1300 1400 1500 1600 1700 1800 1900 2000 
Dilution Rate Input Response 
0.1 
& 
a 0.05 
“7200 1300 1400 1500 1600 1700 1800 1900 2000 
Airflow Rate Input Response 
a Aa. PL Pd de he 
= o 
1200 1300 1400 1500 1600 1700 1800 1900 2000 
Time (hr) 


Figure 1. Closed-loop response of substrate, dissolved oxygen, dilution rate and airflow rate 


5.1. Stability analysis of subspace based DAMPC 

In this section, the stability of the DAMPC is investigated. The design methodology is based on 
receding horizon control, which is also used in a standard model predictive control (MPC) system design. 
However, the main differentiating part between the DAMPC and a commonly used MPC is the starting point 
in the design. The starting point for the usual design of MPC is the availability of a dynamic model that 
describes the system to be controlled. With given sufficient accuracy of the model, long prediction horizon 
and control horizon, the closed-loop stability of MPC is obtained. In contrast, the starting point for the design 
of subspace-based DAMPC is the availability of excitation test data for the system to be controlled. Here, 
the question on the closed-loop stability of DAMPC is associated with the excitation test data for the system, 
as well as the prediction horizon and control horizon. 
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Several simulations of DAMPC showed that large prediction horizons may also resulted in larger 
errors of predictions since the actual system is nonlinear and the prediction based on subspace matrices is 
linear. Also, the longer the prediction horizon, of the DAMPC would approach as the classical MPC. 
The focus of the investigation on stability of the closed-loop system was on the quantity and quality of 
the excitation test data. Particularly, simulation studies are carried out to investigate the effect of both 
the window size of data length and prediction horizon length and the closed-loop performance criterion for 
various numerical values are shown. The best results obtained were n = 400 and H, = 35 for both window 
size length and prediction horizon, respectively. The relationship of closed-loop stability to the predictor 
estimation algorithm is proposed. It can be seen, if the prediction converges to the true value, the overall 
closed-loop system will be equivalent to the standard MPC. From the convergence analysis of this algorithm, 
if a long enough prediction horizon for estimation is used, then the multi-step ahead estimated predictions 
will converge to the true multi-step ahead predictions. This can be proven by looking through the linear 
equation used to calculate the estimated output prediction presented again as: 


Fe = LwWp + Luus £ yf + pE (23) 


where: yf : The natural response, or the effect of the initial conditions is considered to be deterministic since 
it resulted from the multiplication by the past inputs and outputs wp. 


Vee The forced response which is clearly resulted from the effect of the input over the prediction 
horizon which is in fact the impulse response coefficients due to the future input increments. 


c 
aa a a ce Lip 24) 
CAi-1 
D 0 0 Ux 
ms CB D Uk+1 
CAÎ?B CA'3B iD ties 
=> 9f = Du, + 524 CAT Buy, & Lytle (25) 


Due to the nature of this algorithm, it will converge globally, which an important property with 
respects to the typical adaptive algorithms is found in the literature. Here we can conclude that if the future 
inputs remain bounded we can assure the convergence and the stability of the implicit predictive model 
which can only be proven by investigating the values of the impulse coefficients incorporated into prediction 
gain matrices over the prediction horizon. Therefore, by plotting Lọ for every manipulated variable, 
the convergence of the predictive model to the true value over the prediction horizon can be illustrated. 
Figure 2 shows the values of L, for both the future dilution rate manipulated variable and airflow rate 
manipulated variable is bounded over the prediction horizon. Therefore, it ensures the convergence of t 
he estimated output prediction to the true value accordingly. 


a) Lu Prediction Gain Matrix for Dilution Rate b) Lu Prediction Gain Matrix for Airflow Rate 
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Figure 2. The values of L, for prediction gain coefficients of, (a) Future dilution rate, (b) Future airflow rate 
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5.2. Robustness tests for closed-loop DAMPC 

In order to access the closed-loop performance of the applied algorithm, a rigorous simulation 
analysis under different operating conditions and kinetic process parameters were considered and discussed 
in this section. The constant influent rates are considered while taking into accounted the kinetic model 
variations and different operation points. Two different performance criterions were calculated for each 
output (C: dissolved oxygen concentration, and S: substrate concentration) which is the integral square error 
(ISE) and the integral absolute error (LIAE) for the purpose of assessment. In Table 4 listed the ISE and IAE 
of the ASP outputs under the applied DAMPC with two different cases of operating conditions: 

- Case 1: The operating condition in this case is Sin = 200 mgl! and Cin = 0.5 mgl! with 
[IX S C X,]’=[215 55 6 400]" and the varied kinetic parameters for the nonlinear term, 
Hmax, Ks, and Ke. 

- Case 2: The operating condition in this case is Sj, =200mgl™* and Cin =0.5mgl™? with, 
[IX S_C XJ = [217.7896 41.2348 6.1146 435.5791]" and the same varied kinetic 
parameters for the nonlinear term are also considered, Umax, Ks, and Ke. 

Results in Table 4 show the applicability of the controller algorithm to adapt both varied kinetic 
parameters and operating conditions. Therefore, system outputs will be regulated and the input reference 
trajectories will be tracked, irregardless of the conditions imposed on the real plant. Furthermore, no explicit 
model is formulated for the plant and any prior information about initial operating conditions is not needed, 
neither does the kinetic parameters. 


Table 4. Performance test numerical values 








Kinetic ISE IAE Kinetic ISE IAE 

parameters S C S C parameters S C S C 
CASE 1 Mmax = 0.15 132.7 1.133 107.9 16.74 CASE2 © pmax =0.15 150.4 1.262 117.3 17.3 

K, = 100 mgt! K, = 100 mg!" 

K. = 2 mgl K. = 2 mgl" 

Umax = 0.21 680.1 18.64 332.8 59.93 Mmax = 0.21 462.2 13.57 246.2 46.2 

K, = 100 mgt! K, = 100 mg!" 

K.=3 mel! K.=3 mel! 

Umax = 0.21 111.8 1.439 113.4 18.87 Umax = 0.21 114.5 1.514 115.1 19.2 

K, = 150 mel! K, = 150 mel! 

K.=2 mel! K,=2 mel! 





6. CONCLUSION 

In this paper, the design of the adaptive data-driven predictive controller for biological ASP is 
presented. The integration between the robustness and efficiency of subspace-based predictor estimation 
principal (where only geometrical tools are deployed without involving into nonlinear optimization); and the 
power of advanced unconstrained MPC design methodology (where MIMO systems can be controlled 
successfully) contributed to a robust-stable DAMPC for nonlinear ASP. From the control perspective, this 
work can be seen as an important investigation in the field ‘model-free’ control design since it has shown that 
the methodology can cope with the same issues of the ‘model-based’ control methods. This includes the ease 
tuning of MPC parameters particularly for multivariable systems. Although, only unconstrained case is 
considered in this paper, we anticipate that the robustness of the constrained case is naturally achieved. 
In conclusion, the results obtained from the implementation of the presented control methodology on 
nonlinear ASP, was encouraging, and it is proposed that these algorithms show much promise for further 
development and analysis in terms of storage requirements, calculation time, and analytical stability proofs to 
gain wide acceptance in industrial scale applications. 
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